{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>\n",
    "<img src=\"../../img/ods_stickers.jpg\">\n",
    "## Открытый курс по машинному обучению\n",
    "<center>Автор материала: Елена Девятайкина (@elenadevyataykina)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Оценка плотности, используя Gaussian Mixture Models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В данном курсе были рассмаотрены следующие алгоритмы кластеризации: \n",
    "- K-means\n",
    "\n",
    "- Affinity Propagation\n",
    "\n",
    "- Спектральная кластеризация\n",
    "\n",
    "- Агломеративная кластеризация\n",
    "\n",
    "Здесь будет рассмотрен EM-алгоритм, реализованный в модуле `sklearn.mixture.GMM`.\n",
    "\n",
    "**Что это такое?**\n",
    "\n",
    "Это алгоритм, используемый для нахождения оценок максимального правдоподобия параметров вероятностных моделей, в случае, когда модель зависит от некоторых скрытых переменных.\n",
    "\n",
    "**Как он работает?**\n",
    "\n",
    "Оценивая максимальное правдоподобие, EM-алгоритм создает модель, которая назначает метки класса точкам данных.\n",
    "\n",
    "EM-алгоритм начинает с того, что пытается сделать вывод на основании параметров модели.\n",
    "\n",
    "Затем следует итерационный трехшаговый процесс:\n",
    "\n",
    "1. E-шаг: На этом шаге на основании параметров модели вычисляются вероятность принадлежности каждой точки данных к кластеру.\n",
    "2. М-шаг: Обновляет параметры модели в соответствии с кластерным распределением, проведенным на шаге E.\n",
    "3. Предыдущие два шага повторяются до тех пор, пока параметры модели и кластерное распределение не уравняются.\n",
    "\n",
    "Часто данный алгоритм используют для разделения смеси гаусиан (GMM). То есть в результате GMM-EM получаем полноценную смесь Гауссианов полученную в локальном максимуме правдоподобия.\n",
    "\n",
    "Далее рассмаотрим GMM более подробно.\n",
    "\n",
    "\n",
    "**GMM** - один из алгоритмов кластеризации наряду с иерархической кластеризацией (кластеризация на основе связности), k-means (кластеризация на основе центров), DBSCAN (кластеризация на основе плотности). В основе данного алгоритма лежат распределения. \n",
    "\n",
    "**Гауссова смесь распределений (GMM — Gaussian mixture models)** - статистическая модель для представления нормально распределенных подвыборок внутри общей выборки. Гауссова смесь распределений параметризируется двумя типами значений — смесь весов компонентов и средних компонентов или ковариаций (для многомерного случая) и др. Если количество компонентов известно, техника, чаще всего используемая для оценки параметров смеси распределений — ЕМ-алгоритм.\n",
    "\n",
    "Метод GMM непосредственно вытекает из теоремы, гласящей, что любая функция плотности вероятности может быть представлена как взвешенная сумма нормальных распределений\n",
    "\n",
    "$$\n",
    "p = \\sum_{j=1}^k w_i \\varphi (x; \\theta_j),$$\n",
    "\n",
    "$$\n",
    "\\sum_{j=1}^k w_j = 1,$$\n",
    "где $\\varphi (x; \\theta_j)$ - функция распределения многомерного аргумента $x$ с параметрами $\\theta_j,$ \n",
    "\n",
    "$$\\varphi (x; \\theta_j) = p(x | \\mu_j, R_j) = \\frac{1}{{2\\pi}^{n/2}{|R_j|}^{1/2}}e^{-\\frac{1}{2}{(x-\\mu_j)}^TR_j^{-1}(x-\\mu_j)},  x \\in \\mathbb R^n$$\n",
    "\n",
    "$w_j$ - ее вес, $k$ - количество компонент в смеси. Здесь $n$ - размерность пространства признаков, $\\mu_j \\in \\mathbb R^n$ - вектор математического ожидания $j$-й компоненты смеси, $R_j \\in \\mathbb R^{n*n}$ - ковариационная матрица.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "\n",
    "warnings.simplefilter('ignore')\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy import stats\n",
    "from sklearn.mixture import GMM\n",
    "from sklearn.neighbors import KernelDensity\n",
    "\n",
    "import seaborn as sns; sns.set()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Введение в GMM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Создадим одномерный набор данных с азаднным распределением и построим гистограмму"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(2)\n",
    "x = np.concatenate([np.random.normal(0, 2, 2000),\n",
    "                    np.random.normal(5, 5, 2000),\n",
    "                    np.random.normal(3, 0.5, 600)])\n",
    "plt.hist(x, 80, normed=True)\n",
    "plt.xlim(-10, 20);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "С помощью GMM отобразим плотность распределения"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = x[:, np.newaxis]\n",
    "clf = GMM(4, n_iter=500, random_state=3).fit(X)\n",
    "xpdf = np.linspace(-10, 20, 1000)\n",
    "density = np.exp(clf.score(xpdf[:, np.newaxis]))\n",
    "\n",
    "plt.hist(x, 80, normed=True, alpha=0.5)\n",
    "plt.plot(xpdf, density, '-r')\n",
    "plt.xlim(-10, 20);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посмотрим на полученные атрибуты GMM.\n",
    "\n",
    "Среднее значение каждой компоненты:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf.means_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ковариация каждой компоненты"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf.covars_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Вес каждой компоненты"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf.weights_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Так как в GMM задано 4 компоненты, построим все плотности распределения"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(x, 80, normed=True, alpha=0.3)\n",
    "plt.plot(xpdf, density, '-r')\n",
    "\n",
    "for i in range(clf.n_components):\n",
    "    pdf = clf.weights_[i] * stats.norm(clf.means_[i, 0],\n",
    "                                       np.sqrt(clf.covars_[i, 0])).pdf(xpdf)\n",
    "    plt.fill(xpdf, pdf, facecolor='gray',\n",
    "             edgecolor='none', alpha=0.3)\n",
    "plt.xlim(-10, 20);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Для построенной модели можно использовать одно из нескольких средних, чтобы оценить, насколько хорошо он подходит данным.\n",
    "Обычно для этого используется **Информационный критерий Акаике (AIC)** или **Байесовский информационный критерий (BIC)**.\n",
    "\n",
    "**Информационный критерий Акаике** - критерий, применяющийся исключительно для выбора из нескольких статистических моделей.\n",
    "\n",
    "$$aic = 2k - 2l,$$\n",
    "\n",
    "где $k$ - число параметров в статистической модели, $l$ - значение логарифмической функции правдоподобия построенной модели.\n",
    "Критерий не только вознаграждает за качество приближения, но и штрафует за использование излишнего количества параметров модели. Считается, что наилучшей будет модель с наименьшим значением критерия aic.\n",
    "\n",
    "**Байесовский информационный критерий** - критерий, разработанный исходя из баесовского подхода, модификация aic.\n",
    "\n",
    "$$bic = k ln(n) - 2l$$\n",
    "\n",
    "Данный критерий налагает больший штраф на увеличение количества параметров по сравнению с aic."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посмотрим на значение aic, bic на заданной выборке"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(clf.bic(X))\n",
    "print(clf.aic(X))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Давайте рассмотрим их как функцию числа гауссиан:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_estimators = np.arange(1, 10)\n",
    "clfs = [GMM(n, n_iter=1000).fit(X) for n in n_estimators]\n",
    "bics = [clf.bic(X) for clf in clfs]\n",
    "aics = [clf.aic(X) for clf in clfs]\n",
    "\n",
    "plt.plot(n_estimators, bics, label='BIC')\n",
    "plt.plot(n_estimators, aics, label='AIC')\n",
    "plt.xlabel('Number of components')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Судя по графику можно сказать, что для aic предпочтительнее 4 компоненты, а для bic - 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Так как GMM является Generative Model(вероятностная модель для генерации набора данных). С помощью даннной модели можно обнаружить выбросы (Рассчитывается вероятность каждой точки под моделью и чем ниже $p(x)$, тем вероятнее, что данная точка является выбросом. Есть подробная статья, связанная с этой темой - https://arxiv.org/pdf/1706.02690.pdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(0)\n",
    "\n",
    "# Добавим к выборке 20 выбросов\n",
    "true_outliers = np.sort(np.random.randint(0, len(x), 20))\n",
    "y = x.copy()\n",
    "y[true_outliers] += 50 * np.random.randn(20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf = GMM(4, n_iter=500, random_state=0).fit(y[:, np.newaxis])\n",
    "xpdf = np.linspace(-10, 20, 1000)\n",
    "density_noise = np.exp(clf.score(xpdf[:, np.newaxis]))\n",
    "\n",
    "plt.hist(y, 80, normed=True, alpha=0.5)\n",
    "plt.plot(xpdf, density_noise, '-r')\n",
    "plt.xlim(-15, 30);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Далее оценим логарифмическую вероятность каждой точки под моделью"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "log_likelihood = clf.score_samples(y[:, np.newaxis])[0]\n",
    "plt.plot(y, log_likelihood, '.k');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "detected_outliers = np.where(log_likelihood < -9)[0]\n",
    "\n",
    "print(\"Все выбросы:\")\n",
    "print(true_outliers)\n",
    "print(\"\\nОбнаруженные выбросы:\")\n",
    "print(detected_outliers)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Как можно заметить, алгоритм пропустил несколько точек (те, которые находятся в середине распределения).\n",
    "\n",
    "Вот пропущенные выбросы:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "set(true_outliers) - set(detected_outliers)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "И вот значения, которые были неверно определены как выбросы"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "set(detected_outliers) - set(true_outliers)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Необходимо отметить, что здесь используются одномерные данные, но GMM также может делать обобщения на несколько измерений."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Другие оценки плотности"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Важным средством оценки плотности является оценка плотности ядра, реализованная в `sklearn.neighbors.KernelDensity`. Можно это рассматривать как обобщение GMM. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kde = KernelDensity(0.15).fit(x[:, None])\n",
    "density_kde = np.exp(kde.score_samples(xpdf[:, None]))\n",
    "\n",
    "plt.hist(x, 80, normed=True, alpha=0.5)\n",
    "plt.plot(xpdf, density, '-b', label='GMM')\n",
    "plt.plot(xpdf, density_kde, '-r', label='KDE')\n",
    "plt.xlim(-10, 20)\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Все оценки плотности можно рассматривать как генерирующие модели данных: то есть это позволяет понять, как создавать больше данных, которые соответствуют модели. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
